Tanzania, August, 2019

Introduction to leaflet

Intro to leaflet

What is leaflet?

  • open source JavaScript library used to build web mapping applications
  • Probably the world's most popular web mapping system

Why leaflet?

  • Great way to explore datasets with geographic component
  • Communicate results with interactive maps using only a few lines of code
  • Integrates nicely with Rshiny producing powerful extra-interactive maps

Intro to leaflet

Getting started with our map

leaflet() initialises a leaflet object m. Pipe this into addCircles() where we refer to the map_data dataframe and indicate which columns correspond to latitude and longitude.

m <- leaflet() %>%
  addCircles(data = map_data, lng = ~x, lat = ~y)
m

Tiles

Tiles give our data geographic context. addTiles() uses OpenStreetMap tiles, the default for leaflet

m %>% addTiles()

Using addProviderTiles()

Allows us to use third party tiles (e.g. Terrain tiles from Stamen Design).

# m created using code from previous slides
m <- leaflet() %>%
  addCircles(data=map_data, lng=~x, lat=~y)

Code options

addProviderTiles(map = m, provider = "Stamen.Terrain")

addProviderTiles(map = m, provider = providers$Stamen.Terrain)

m %>% addProviderTiles(provider = providers$Stamen.Terrain)

Your turn 2.2a

Try some options for provider tiles.

Hint: Use tab completion with addProviderTiles() to access options for arguments and available tiles.

Hint: providers$

Using addProviderTiles()

m %>% addProviderTiles(providers$Esri.NatGeoWorldMap)

Using addProviderTiles()

m %>% addProviderTiles(providers$Wikimedia)

Using addProviderTiles()

m %>% addProviderTiles(providers$Esri.WorldTopoMap)

Using addProviderTiles()

m %>% addProviderTiles(providers$Esri.WorldImagery)

Using addProviderTiles()

m %>% addProviderTiles(providers$CartoDB.Positron)

More info on addProviderTiles()

Popups

An interactive way for the user to get more info on data points

Easier to first store what we want to display as an object popupInfo

popupInfo <- paste("Date: ", map_data$date, "<br>",
                   "Species: ", map_data$species, "<br>",
                   "Age: ", map_data$age, "<br>",
                   sep = " ")

m <- leaflet() %>% addProviderTiles(providers$CartoDB.Positron)

Then, refer to popupInfo when using addCircles(), addCircleMarkers() etc.

m %>% addCircles(data = map_data, lng = ~x, lat = ~y, popup = popupInfo)

Popups

An interactive way for the user to get more info on data points

m %>% addCircles(data = map_data, lng = ~x, lat = ~y, popup = popupInfo)

Popups and labels

Labels show up when hovered over, popups when the user clicks.

m %>% addCircles(data = map_data, lng = ~x, lat = ~y, popup = popupInfo,
                 label = ~species)

Circle markers

Circle markers (addCircleMarkers()) are similar to circles but radius in onscreen pixels is constant regardless of zoom level.

m %>% addCircleMarkers(data=map_data, lng=~x, lat=~y, popup = popupInfo, radius = 5,
                       stroke = T, weight = 3, color = "grey", opacity = 0.8,
                       fillColor = ifelse(map_data$species == "human",
                                          "firebrick", "dodgerblue"), fillOpacity = 0.7)

Circle markers

With, stroke = F, weight, color, and opacity are irrelevant. Can also use radius = to size circle markers by species.

m %>% addCircleMarkers(data = map_data, lng = ~x, lat = ~y, popup = popupInfo, stroke = F,
                       radius = ifelse(map_data$species == "human", 8, 4),
                       fillColor = ifelse(map_data$species == "human", "firebrick", "dodgerblue"),
                       fillOpacity = 0.8)

addMarkers

On a map with many markers, it might make sense to cluster them using clusterOptions = markerClusterOptions() shown here with the default options for markerClusterOptions().

m %>% addMarkers(data = map_data, lng = ~x, lat = ~y, label = ~species, popup = popupInfo)

addMarkers

On a map with many markers, it might make sense to cluster them using clusterOptions = markerClusterOptions() shown here with the default options for markerClusterOptions().

m %>% addMarkers(data = map_data, lng = ~x, lat = ~y, label = ~species, popup = popupInfo,
                 clusterOptions = markerClusterOptions())

Colour palettes in leaflet

Discrete variable - using colorFactor()

leaflet::colourFactor() creates a function (we're going to call it myPal) that is used as an argument when we use addCircles(), addCircleMarkers()etc.

# Choose 3 colours for our 3 groups  
colour_pal <- c("gold", "forestgreen", "cornflowerblue")
myPal <- colorFactor(palette = colour_pal,
                     domain = c("Human", "Domestic", "Wildlife"))

Create a fresh version of m before we add coloured circles

m <- leaflet() %>% 
  addProviderTiles(provider = providers$CartoDB.Positron)

Discrete variable - using colorFactor()

myPal() takes species type variable as argument

m <- m %>% addCircles(data = map_data, lng = ~x, lat = ~y, 
                      color = myPal(map_data$species_type),
                      opacity = 0.9, popup = popupInfo)
m

Adding a legend

m <-  m %>% addLegend(position = "topright", title = "Species<br>type", pal = myPal,
                      values = map_data$species_type, opacity = 0.9)
m

Your turn 2.2b

Adapt the code from previous slides to produce a map with:
1) tiles showing topography
2) circles coloured by species (instead of species type)
3) a colour legend
4) a scale bar showing distance in kilometers (hint: ?addScaleBar).

Previous code:

colour_pal <- c("gold", "forestgreen", "cornflowerblue")
myPal <- colorFactor(colour_pal, domain = c("Human", "Domestic", "Wildlife"))
m <- leaflet() %>%
  addProviderTiles(provider = providers$CartoDB.Positron) %>%
  addCircles(data = map_data, lng = ~x, lat = ~y,
             color = myPal(map_data$species_type), opacity = 0.9) %>%
  addLegend(position = "topright", title = "Species<br>type", pal = myPal,
            values = map_data$species_type, opacity = 0.9)
m

Solution 2.2b

Part 1) tiles showing topography

m <- leaflet() %>%
  addProviderTiles(provider = providers$Esri.WorldPhysical)
m %>% setView(lat = -6.4333, lng = 38.9, zoom = 5)

Solution 2.2b

Part 2) circles coloured by species (instead of species type)

colour_pal <- c("gold", "forestgreen", "cornflowerblue", "red", "magenta")
myPal <- colorFactor(palette = colour_pal, domain = unique(map_data$species))
m <- m %>% addCircles(data = map_data, lng = ~x, lat = ~y,
                      color = myPal(map_data$species), opacity = 0.9)
m

Solution 2.2b

Part 2) circles coloured by species (instead of species type)

myPal <- colorFactor(palette = "Set2", domain = unique(map_data$species))
m <- m %>% addCircles(data = map_data, lng = ~x, lat = ~y,
                      color = myPal(map_data$species), opacity = 0.9)
m

Solution 2.2b

Part 3) a colour legend

m <- m %>% addLegend(position = "topright", title = "Species", pal = myPal,
                     values = map_data$species, opacity = 0.9)
m

Solution 2.2b

Part 4) a scale bar showing distance in kilometers

m <- m %>% addScaleBar(position = "bottomleft",
                       options = scaleBarOptions(maxWidth = 200, # default = 100
                                                 imperial = FALSE))
m

Continuous variable - colour by date

First, need to convert date variable to a numeric value

head(map_data$date, 5)
## [1] "2014-01-01" "2014-01-01" "2014-01-02" "2014-01-03" "2014-01-04"
is.factor(map_data$date)
## [1] FALSE
# Use 2 functions from lubridate for conversion
map_data$date_decimal <- lubridate::decimal_date(lubridate::ymd(map_data$date))

Continuous variable - colour by date

date_decimal is now a numeric variable

# Use 2 functions from lubridate for conversion
map_data$date_decimal <- lubridate::decimal_date(lubridate::ymd(map_data$date))

head(map_data$date_decimal)
## [1] 2014.000 2014.000 2014.003 2014.005 2014.008 2014.014
is.numeric(map_data$date_decimal)
## [1] TRUE

Continuous variable - colour by date

Provide some colours and values of date to colorNumeric().

Either provide colours (palette = c("darkgreen", "green", "gold""))) or refer to a palette (palette = "Spectral").

dateRange <- c(2014, 2015)
myPal <- colorNumeric(palette = "YlOrRd", domain = dateRange, reverse = T)

# new leaflet widget
m <- leaflet() %>% addProviderTiles(providers$CartoDB.Positron)

Continuous variable - colour by date

Add circles and a legend referring to myPal

m %>% addCircles(data = map_data, lng = ~x, lat = ~y,
                 color = myPal(map_data$date_decimal), opacity = 0.9) %>% 
  addLegend(position = "topright", title = "Date",
            pal = myPal, values = dateRange, opacity = 0.9,
            labFormat = labelFormat(big.mark = "")) # removes a comma from 2,014

Your turn 2.2c

Adapt the code from previous slides to produce a map with:
1) Circles coloured by population density using the Spectral palette
2) and sized by population density (higher density = larger circles)
3) with colour legend and scale bar

If you get the solution, try colorBin() instead of colourNumeric()

Previous code to adapt:

dateRange <- c(2014, 2015)
myPal <- colorNumeric(palette = "YlOrRd", domain = dateRange)
m <- leaflet() %>%
  addProviderTiles(provider = providers$CartoDB.Positron) %>%
  addCircles(data = map_data, lng = ~x, lat = ~y,
             color = myPal(map_data$date_decimal), opacity = 0.9) %>% 
  addLegend(position = "topright", title = "Date",
            pal = myPal, values = dateRange, opacity = 0.9,
            labFormat = labelFormat(big.mark = ""))
m

Solution 2.2c

Colouring by density:

densityRange <- c(0.1, 12000)
myPal <- colorNumeric(palette = "Spectral", domain = densityRange)

m <- leaflet() %>%
  addProviderTiles(provider = providers$CartoDB.Positron) %>%
  addCircleMarkers(data = map_data, lng = ~x, lat = ~y,  radius = 3,
             stroke = F, fillColor = myPal(map_data$density), fillOpacity = 0.9) %>% 
  addLegend(position = "topright", title = "People<br>per km<sup>2</sup>",
            pal = myPal, values = densityRange, opacity = 0.9) %>% 
  addScaleBar(position = "bottomleft")

Solution 2.2c

Colouring by density:

m

Solution 2.2c

and sizing by density, we add radius = ~log(density):

densityRange <- c(0.1, 12000)
myPal <- colorNumeric(palette = "Spectral", domain = densityRange)

m <- leaflet() %>%
  addProviderTiles(provider = providers$CartoDB.Positron) %>%
  addCircleMarkers(data = map_data, lng = ~x, lat = ~y, radius = ~log(density),
             stroke = F, fillColor = myPal(map_data$density), fillOpacity = 0.9) %>% 
  addLegend(position = "topright", title = "People<br>per km<sup>2</sup>",
            pal = myPal, values = densityRange, opacity = 0.9) %>% 
  addScaleBar(position = "bottomleft")

Solution 2.2c

and sizing by density, we add radius = ~log(density):

m

Solution 2.2c with colorBin()

With colorBin() we provide a vector of numbers

When we build up the leaflet widget, we just remove the log() functions as we no longer need to transform

myPal <- colorBin(palette = "Spectral", domain = map_data$density,
                  bins = c(0, 100, 250, 500, 750, 1000, 5000, Inf))

m <- leaflet() %>%
  addProviderTiles(provider = providers$CartoDB.Positron) %>%
  addCircleMarkers(data = map_data, lng = ~x, lat = ~y, radius = ~log(density),
             stroke = F, fillColor = myPal(map_data$density), fillOpacity = 0.9) %>% 
  addLegend(position = "topright", title = "Human<br>density",
            pal = myPal, values = map_data$density, opacity = 0.9) %>% 
  addScaleBar(position = "bottomleft")

Solution 2.2c

m

colourQuantile()

myPal <- colorQuantile(palette = "Spectral", domain = map_data$density, n = 5)

leaflet() %>%
  addProviderTiles(provider = providers$CartoDB.Positron) %>%
  addCircleMarkers(data = map_data, lng = ~x, lat = ~y, radius = ~log(density),
             stroke = F, fillColor = myPal(map_data$density), fillOpacity = 0.9) %>%
  addLegend(position = "topright", title = "Human<br>density",
            pal = myPal, values = map_data$density, opacity = 0.9,
            labFormat = labelFormat(big.mark = ""))

Using leaflet and shape file data

Shape file data

Load shape data and look at the structure

region_shp <- rgdal::readOGR(dsn = "data/TZ_Region_2012_density")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/wharvey/Documents/git/ShinyCourse/Day 2/slides/data/TZ_Region_2012_density", layer: "TZ_Region_2012.density"
## with 30 features
## It has 4 fields
head(region_shp)
## class       : SpatialPolygonsDataFrame 
## features    : 6 
## extent      : 31.03802, 39.85569, -6.976504, -1.027046  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +no_defs 
## variables   : 4
## names       :   Region_Nam,              Loc_ID, Loc_type,    density 
## min values  : Kusini Pemba, Region-Kusini Pemba,   Region,   94.02598 
## max values  :        Tanga,        Region-Tanga,   Region, 1002.31128

Plotting polygons

Use addPolygons() to add to leaflet map. label = determines what will be shown when cursor hovers above area.

m <- leaflet() %>% addProviderTiles(provider = providers$Esri.WorldShadedRelief)
m %>% addPolygons(data = region_shp, label = region_shp$Region_Nam,
                  color = "black", weight = 1, fillOpacity = 0.3, opacity = 0.6)

Plotting polygons

Use highlightOptions() to highlight shapes when cursor hovers above. bringToFront = TRUE avoids potential problems with order in which shapes are coded in shape file.

m %>% addPolygons(data = region_shp, label = region_shp$Region_Nam, fillColor = "gold",
                  color = "black", weight = 1, fillOpacity = 0.3, opacity = 0.6,
                  highlightOptions = highlightOptions(color = "white", weight = 3,
                                                      bringToFront = TRUE, opacity = 1))

Plotting polygons

Possible to use addPolylines() to add outlines only to leaflet map.

m %>% addPolylines(data = region_shp, color = "red",  weight = 2, opacity = 1)

Plotting a choropleth

Choropleth: Map with areas shaded according to some numeric variable
(similar to heatmap but uses geographic boundaries)

myPal <- colorNumeric(palette = "YlOrRd", domain = region_shp$density)
m %>% addPolygons(data = region_shp, label = region_shp$Region_Nam,
                  fillColor = myPal(region_shp$density), fillOpacity = 0.7,
                  color = "white", opacity = 1, weight = 2)

Plotting a choropleth

Like before, we can log density to give a more informative map

myPal <- colorNumeric(palette = "YlOrRd", domain = log(region_shp$density))
m %>% addPolygons(data = region_shp, label = region_shp$Region_Nam,
                  fillColor = myPal(log(region_shp$density)), fillOpacity = 0.7,
                  color = "white", opacity = 1, weight = 2)

Your turn 2.2d

Adjust the choropleth code below to:

  1. use colorQuantile() instead of colorNumeric()
  2. add circles coloured by date for human cases
m <- leaflet() %>% addProviderTiles(provider = providers$Esri.WorldShadedRelief)

myPal <- colorNumeric(palette = "YlOrRd", domain = log(region_shp$density))
m %>% addPolygons(data = region_shp, label = region_shp$density,
                  fillColor = myPal(log(region_shp$density)), fillOpacity = 0.7,
                  color = "white", opacity = 1, weight = 2)

Solution 2.2d

Subset using map_data[map_data$species == "human",] or alternatively dplyr::select()

With colorQuantile(), we don't need to log density

Create a second colour palette to colour cases by date using colorNumeric().

human_data <- map_data[map_data$species == "human",]
myPal <- colorQuantile(palette = "YlOrRd", domain = region_shp$density)
myPal2 <- colorNumeric(palette = "viridis", domain = human_data$date_decimal)
m %>% addPolygons(data = region_shp, label = region_shp$density,
                  fillColor = myPal(region_shp$density), fillOpacity = 0.7,
                  color = "white", opacity = 1, weight = 2) %>% 
  addCircleMarkers(data = human_data, lat = ~y, lng = ~x, stroke = T,
                   fillColor = myPal2(human_data$date_decimal), fillOpacity = 0.7,
                   radius = 4, weight = 2, color = "black", opacity = 0.8)

Solution 2.2d

human_data <- map_data[map_data$species == "human",]
myPal <- colorQuantile(palette = "YlOrRd", domain = region_shp$density)
myPal2 <- colorNumeric(palette = "viridis", domain = human_data$date_decimal)
m %>% addPolygons(data = region_shp, label = region_shp$density,
                  fillColor = myPal(region_shp$density), fillOpacity = 0.7,
                  color = "white", opacity = 1, weight = 2) %>% 
  addCircleMarkers(data = human_data, lat = ~y, lng = ~x,
                   fillColor = myPal2(human_data$date_decimal), fillOpacity = 0.7, radius = 4,
                   stroke = T, weight = 2, color = "black", opacity = 0.8)

Add labels with some html formatting

labels <- sprintf("<strong>%s</strong><br/>%g people / km<sup>2</sup>",
                   region_shp$Region_Nam, round(region_shp$density, 0)) %>%
  lapply(htmltools::HTML)

m %>% addPolygons(data = region_shp, fillColor = myPal(region_shp$density), label = labels,
                  weight = 2, opacity = 1, color = "white", fillOpacity = 0.7,
                  highlightOptions = highlightOptions(color = "black", weight = 3,
                                                      bringToFront = TRUE, opacity = 1))

Using leaflet and environmental data in raster format

Loading raster data

Load data using raster() from the raster package

density <- raster::raster("data/HumanDensity/HumanPopulation.grd")
density
## class       : RasterLayer 
## dimensions  : 259, 267, 69153  (nrow, ncol, ncell)
## resolution  : 0.041665, 0.041665  (x, y)
## extent      : 29.32649, 40.45105, -11.77259, -0.9813543  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/wharvey/Documents/git/ShinyCourse/Day 2/slides/data/HumanDensity/HumanPopulation.grd 
## names       : TZA_popmap10adj_v2b 
## values      : 0, 29089.56  (min, max)

Plotting raster data

We can use the base plotting function plot() on the object we've created

plot(density)

Plotting raster data

Logging density will provide more texture to map

plot(log(density))

Plotting raster data with leaflet

We use colorNumeric() again, as we did to colour by date earlier

colour_pal <- c("#FFFFCC", "#41B6C4", "#0C2C84")
raster::values(density) %>% range(na.rm = T)
## [1]     0.00 29089.56
myPal <- colorNumeric(palette = colour_pal, domain = c(0, 30000), na.color = "red")

m <- leaflet() %>%
  addRasterImage(x = density, colors = myPal, opacity = 0.8)

Plotting raster data with leaflet

m

Plotting raster data with leaflet

Let's log density as we did before

raster::values(density) %>% range(na.rm = T) %>% log
## [1]     -Inf 10.27813
myPal <- colorNumeric(palette = colour_pal, domain = c(-1.5, 11),
                      na.color = "transparent")

We have to exclude some data, as there are 0 density polygons and log(0) gives -infinity

sum(raster::values(density) > exp(-1.5), na.rm = T) /
  sum(raster::values(density) >= 0, na.rm = T)
## [1] 0.9975207

Plotting raster data with leaflet

Plot the logged raster data with some tiles and a legend for the colour scheme.

m <- leaflet() %>%
  addProviderTiles(provider = providers$CartoDB.Positron) %>%
  addRasterImage(x = log(density), colors = myPal, opacity = 0.8) %>%
  addLegend(pal = myPal, values = c(-1.5, 11), opacity = 0.8,
            title = "Human<br>population<br>density")
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA

Plotting raster data with leaflet

Plotting raster data with leaflet

Instead of using addProviderTiles(provider = providers$CartoDB.Positron)

We order things here so that we plot

  1. The background tiles only (providers$CartoDB.PositronNoLabels)
  2. Then the raster data (addRasterImage())
  3. and finally the labels (providers$CartoDB.PositronOnlyLabels)
m <- leaflet() %>%
  addProviderTiles(provider = providers$CartoDB.PositronNoLabels) %>%
  addRasterImage(x = log(density), colors = myPal, opacity = 0.8) %>%
  addProviderTiles(provider = providers$CartoDB.PositronOnlyLabels) %>%
  addLegend(pal = myPal, values = c(-1.5, 11), opacity = 0.8,
            title = "Human<br>population<br>density")

Plotting raster data with leaflet

Your turn 2.2e

Adapt the code below to add region outlines that highlight when the cursor hovers an area. Use addPolygon() with fillOpacity = 0

colour_pal <- c("#FFFFCC", "#41B6C4", "#0C2C84")
myPal <- colorNumeric(palette = colour_pal, domain = c(-1.5, 11),
                      na.color = "transparent")

m <- leaflet() %>%
  addProviderTiles(provider = providers$CartoDB.PositronNoLabels) %>%
  addRasterImage(x = log(density), colors = myPal, opacity = 0.8) %>%
  addProviderTiles(provider = providers$CartoDB.PositronOnlyLabels) %>%
  addLegend(pal = myPal, values = c(-1.5, 11), opacity = 0.8,
            title = "Human<br>population<br>density")
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA

Solution 2.2e

m <- leaflet() %>%
  addProviderTiles(provider = providers$CartoDB.PositronNoLabels) %>%
  addRasterImage(x = log(density), colors = myPal, opacity = 0.8) %>%
  addProviderTiles(provider = providers$CartoDB.PositronOnlyLabels) %>%
  addLegend(pal = myPal, values = c(-1.5, 11), opacity = 0.8,
            title = "Human<br>population<br>density") %>% 
  addPolygons(data = region_shp, label = region_shp$Region_Nam, fillOpacity = 0,
              color = "black", weight = 1, opacity = 0.6,
              highlightOptions = highlightOptions(color = "white", weight = 3,
                                                  bringToFront = TRUE, opacity = 1))

Solution 2.2e

Exporting leaflet map

Exporting as a widget

saveWidget() from the htmlwidgets package can save the map as a .html file.

You can email this to someone, add to a webpage etc.

There is no need for R to be running in the background.

saveWidget(widget = m, file = "filename.html")

Exporting as an image

mapshot() from the mapview package can save the map as a static image in either a .png, .pdf, or .jpeg file.

m <- m %>% setView(m, lat = -6.4333, lng = 38.9, zoom = 7.5)
mapshot(m, file = "filename.png", width = 7, height = 6)

Exporting as an image

mapshot() from the mapview package can save the map as a static image in either a .png, .pdf, or .jpeg file.

Combine with setView() from the leaflet package to control where map view is centered and level of zoom

m <- m %>% setView(lat = -6.4333, lng = 38.9, zoom = 7.5)
mapshot(m, file = "filename_zoom.png", width = 7, height = 6)